An $n \times n$ matrix $A$ is invertible, if and only if $A$ is row equivalent to $I_n$, the $n \times n$ identity matrix. In this case, the row operations that reduces $A$ to $I_n$, reduces $I_n$ to $A^{-1}$.
You will find the invert of $A$ by applying the Gauss-Jordan elimination to the aumented matrix $ [A|I_n]$.
In order to reuse the code you wrote for your last homeworks, we will split the computation in two steps:
and a upward pass of Gaussian elimination that will reduce your augmented matrix to a row echelon form.
We will provide you with some useful functions and scripts to test and debug your functions.
In [ ]:
function rowSwap!(A, i,j)
# input: A matrix
# i,j the row indexes to be swapped
n = size(A)[1]
# checking that i and j are within the permitted range
(i > n || j > n ) && error("Index out of range in the rowSwap")
# if the index are not the same
if i != j
buffer = A[i,:]
A[i,:] = A[j,:]
A[j,:] = buffer
end
end
Q1.a You will write a function that performs the same Gaussian elimination that you wrote for homework 4 but with one difference. The function takes as input:
In other words, your matrix will take a matrix as a right-hand-side.
Your function will create the augmented system, and perform the reduction to a upper triangular matrix.
The output of the function will be the tuple (U,B1).
U is the upper triangular matrix resulting from the elimination and $B1$, is the second block of the augmented matrix.
To obtain $U$ and $B1$ from your augmented matrix you perform a slicing (i.e. use [:,1:n] and [:,n+1:end]).
You will use your rowSwap function acting on the augmented matrix.
Finally, your function will raise an error if:
Hints:
In [ ]:
function gaussianElimination(A,B)
#input: A squared matrix
# b a vector
#output: U upper triangular matrix
# b1 the resulting vector
# safety checks
(n,m) = size(A)
(n != m) && error("Matrix is not square \n")
(n != size(B)[1]) && error("Dimension mismatch \n")
# create the augmented matrix
M = hcat(A,B)
for i = 1:n
# look for the first non zero entry in the ith column
# find the indices of the non zero elements (up to machine precision)
indexj = find(abs(M[i:end,i]).> eps(1.0))
# if indexj is empty then the matrix is singular and raise error
isempty(indexj) && error("The matrix is singular \n")
# call row swap
rowSwap!(M,i,(i-1)+indexj[1])
# for loop for eliminating unknows
for j = i+1:n
M[j,:] = M[j,:] - (M[j,i]/M[i,i])*M[i,:]
end
end
# checking the last pivot!!
abs(M[n,n]) < eps(1.0) && error("The matrix is singular \n")
# slicing the matrices
U = M[:,1:n]
B1 = M[:,n+1:end]
return (U,B1)
end
Q1.b You will write a function that perform the second part of the Gauss-Jordan method, i.e. the reduction to reduced row echelon form .
The input of your function will be:
You will reduce the augmented matrix $U|B$ using the Gauss-Jordan method. Given that $U$ is already in upper triangular form, you need to reduce it to diagonal form. I.e. you need to eliminate the non-zero elements above the diagonal. In this case the final diagonal form will be an identity matrix.
This transformation to reduced row echelon form is equivalent to a triangular solve with mulitple right-hand-sides. The ouput of the function will be the tuple $(I,X)$. $I$ is the first block of the the augmented matrix, and it should be almost an identity matrix (up to machine precission). $X$ will be the solution to the matrix system $UX= B$. Note that in this case the solution $X$ is a matrix.
Your function needs to have safeguards against a size mismatch (i.e., the sizes of the matrices are not compatible, or your matrix $U$ is not a square matrix).
Hint: if you need to run a foor loop that goes from n-1 to 1, you can use the syntax
for j = n-1:-1:1
In [ ]:
function upwardGaussianElimination(A,B)
#input: A squared matrix
# b a vector
#output: U upper triangular matrix
# b1 the resulting vector
# safety checks
(n,m) = size(A)
(n != m) && error("Matrix is not square \n")
(n != size(B)[1]) && error("Dimension mismatch \n")
# create the augmented matrix
M = hcat(A,B)
for i = n:-1:2
# for loop for eliminating unknows
M[i,:] = M[i,:]./M[i,i]
for j = i-1:-1:1
M[j,:] = M[j,:] - (M[j,i]/M[i,i])*M[i,:]
end
end
M[1,:] = M[1,:]./M[1,1]
# slicing the matrices
I = M[:,1:n]
X = M[:,n+1:end]
return (I,X)
end
You can test that your function is correct by running the following script:
In [ ]:
# size of the Matrix
m = 100
# creating an upper triangular Matrix
U = diagm(m*ones(m,1)[:], 0) + diagm(rand(m-1,1)[:], 1) + diagm(rand(m-2,1)[:], 2)
# creating the rhs
B = rand(size(U))
@time (I,X) = upwardGaussianElimination(U,B)
print("Residual of the backward substitution is ", norm(U*X -B)/norm(B),"\n")
print("Distance to an Identity matrix is ", norm(I - eye(m)), "\n")
The residual should be extremely small (around epsilon machine)
Q1.c You will write a function (very short) that finds the inverse of a Matrix $A$.
The input of your function will be :
Your function will apply the Gaussian Elimination to the augmented matrix $[A|I_n]$, and it will use the Gauss-Jordan method to find the inverse of $A$. Your function needs to check that your matrix is squared.
In [ ]:
function invert(A)
# input: A square matrix
# output: A^{-1}
(U,B) = gaussianElimination(A,eye(size(A)[1]))
(I, Ainv)= upwardGaussianElimination(U,B)
return Ainv
end
You can test your code by running the following script
In [ ]:
# size of the Matrix
m = 100
# creating the Matrix
A = randn(m,m) + m*eye(m)
# compute the inverse
Ainv = invert(A)
# test if Ainv is the inverse
println("Distance to the identity is ", norm(A*Ainv - eye(m)))
The residual should be extremely small (around epsilon machine)
Q2.a You will write a function that performs the Gaussian elimination with partial pivoting. You will as a base the function you wrote for Homework 4, and you will modify it slightly. The function takes as input:
Your function will create the augmented system, and perform the Gaussian elimination.
The output of the function will be the tuple (U,b1).
U is the upper triangular matrix resulting from the elimination and b1, is the resulting vector.
To obtain $U$ and $b1$ from your augmented matrix you perform a slicing (i.e. use [:,1:n]).
You will use your rowSwap function acting on the augmented matrix to perform the partial pivoting.
Finally, your function will raise an error if:
Hints:
In [ ]:
function gaussianEliminationPartialPivoting(A,b)
#input: A squared matrix
# b a vector
#output: U upper triangular matrix
# b1 the resulting vector
# safety checks
(n,m) = size(A)
(n != m) && error("Matrix is not square \n")
(n != size(b)[1]) && error("Dimension mismatch \n")
# create the augmented matrix
M = hcat(A,b)
for i = 1:(n-1)
# look for the first non zero entry in the ith column
# find the indices of the non zero elements (up to machine precision)
indexj = sortperm(M[i:end,i], by=abs)[end] + (i-1)
# if indexj is empty then the matrix is singular and raise error
(abs(M[indexj,i]) < eps(1.0)) && error("The matrix is singular \n")
# call row swap
rowSwap!(M,i,indexj)
# for loop for eliminating unknows
M[i,:] = M[i,:]/M[i,i]
for j = i+1:n
M[j,:] = M[j,:] - M[j,i]*M[i,:]
end
end
abs(M[n,n]) < eps(1.0) && error("The matrix is singular \n")
#normalizing the last row
M[n,:] = M[n,:]/M[n,n]
# slicing the matrices
U = M[:,1:n]
b1 = M[:,n+1:end]
return (U,b1)
end
Q2.b You will use the triangular solver you wrote in the last homework to implement a linear solver.
In [ ]:
# write your triangular solver here
function backwardSubstitution(U,b)
# input: U upper triangular matrix
# b vector
# output: x = U\b
# checks for sizes
(n,m) = size(U)
(n != m) && error("Upper triangular matrix is not square \n")
(n != size(b)[1]) && error("Dimension mismatch \n")
# creating x and running the backward substitution
x = zeros(b)
x[n] = b[n]/U[n,n]
for i = (n-1):-1:1
x[i] = (b[i] - dot(U[i,i+1:end][:],x[i+1:end]))/U[i,i]
end
# returning x
return x
end
In an analogous manner to Homework 4, you will use the Gaussian Elimination algorithm with partial pivoting and your triangular solve to obtain a more stable Linear solver.
In [ ]:
function solveGaussPivoted(A,b)
# input: A square matrix
# b vector
# output: x = A\b
(U,b1) = gaussianEliminationPartialPivoting(A,b)
return backwardSubstitution(U,b1)
end
You will test your function by finding the solution to $Ax =b$ for a mildly ill-conditioned matrix.
In [ ]:
# size of the Matrix
m = 100
# creating the Matrix
A = randn(m,m)
println("The conditioning number of A is " ,cond(A))
# creating the rhs
b = rand(size(A)[1],1)
@time x = solveGaussPivoted(A,b)
print("Residual of the solver is ", norm(A*x -b)/norm(b),"\n")
Q2.c You can use the Gaussian Elimination you wrote in Q1 (without pivoting) to solve a linear system
In [ ]:
function solveGauss(A,b)
# input: A square matrix
# b vector
# output: x = A\b
(U,b1) = gaussianElimination(A,b)
return backwardSubstitution(U,b1)
end
You will compare the pivoted versus the non pivoted solvers. By running the following script.
In [ ]:
# size of the Matrix
m = 100
# creating the Matrix
A = randn(m,m)
println("The conditioning number of A is " ,cond(A))
# creating the rhs
b = rand(size(A)[1],1)
@time xPiv = solveGaussPivoted(A,b)
@time xNonPiv = solveGauss(A,b)
print("Residual of the solver is ", norm(A*xPiv -b)/norm(b),"\n")
print("Residual of the solver is ", norm(A*xNonPiv -b)/norm(b),"\n")
As you can see, each time you execute the script the pivoted and non-pivoted algorithms give you back an answer with different residuals. Why do you see this behavior? How does the pivoting help the stability?
Answer: We can observe that the residual is about two orders of magnitude smaller, this is mainly due to the partial pivoting, which avoids divisions by small pivots.
In [ ]: